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We introduce a set of coupled equations for multilayer water waves that removes the 
ill-posedness of the multilayer Green-Naghdi (MGN) equations in the presence of 
shear. The new well-posed equations are Hamiltonian and in the absence of imposed 
background shear they retain the same travelling wave solutions as MGN. We call 
the new model the Square Root Depth (vZ)) equations, from the modified form of 
their kinetic energy of vertical motion. Our numerical results show how the \AD 
equations model the effects of multilayer wave propagation and interaction, with 
and without shear. 



1. Introduction 

The propagation and interactions of internal gravity waves on the ocean thermo- 
cline may be observed in many areas of strong tidal flow, including the Gibraltar 
Strait and the Luzon Strait. These waves a re strongly nonlinear and may even be 
seen from the Space Shuttle dLiu et aQ 1998). with their crests moving in great arcs 
hundreds of kilometres in length and traversing sea basins thousands of kilome- 
tres across. The M GN model — the multilayer extension of the w ell known Green- 
Naghdi equations dGreen & Naghdill 19761 : IChoi & Camassall 19991) — has been used 



with some success to model the short term behaviour of these waves dJo & Choil 



2002). Nevertheless, MGN and its rigid-lid vers ion, the Choi-Camassa equation 
d 19961. hereafter CC), have both been shown bv lLiska & WendrofJ d 19971) to be ill- 
posed in the presence of background shear. That is, background shear causes the 
linear growth rate of a perturbation to increase without bound as a function of wave 
number. Needless to say, this ill-posedness has made numerical modelling of MGN 
problematic. In particular, ill-posedness prevents convergence of the numerical so- 
lution, since the energy cascades to smaller scales and builds up at the highest 
resolved wave number. Grid refinement only makes the p roblem worse. Regul ariza- 
tion by keeping higher-order expansion terms is possible (Barros & Choill20 09). but 
such methods tend to destroy the Hamiltonian property of the system and thus may 
degrade its travelling wave structure. Nevertheless, if one is to consider the wave 
generation problem, one must consider the effects of both topography and shear. 

Thus, the MGN equations must be modified to make them well-posed. We shall 
require that the new system: 

(a) is both linearly well-posed and Hamiltonian; 

(b) preserves the MGN linear dispersion relation for fluid at rest; 

(c) has the same travelling wave solutions as MGN in the absence of imposed 
background shear. 

The \f~D equation we introduce here satisfies these requirements. The remainder of 
the paper is organised as follows. In Section 2 we derive the \f~D model in the same 
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Euler-Poincare variational framework as for MGN dPercival et aZ.l l2008). In Section 

3 we compare the linear dispersion analysis of \[~D and MGN, and thus show that 
linear ill-posedness has been removed. In Section 4 we show that the \[T) and the 
shear free MGN equations possess the same travelling wave solutions. In Section 
5 we present numerical results that compare the wave propagation and interaction 
properties of the \f~D and GN solutions for a single layer. Finally, in Section 6 we 
present numerical results for the \[T) equations that model the effects of two-layer 
wave propagation and interaction, both with and without a background shear. 



2. The VD Governing Equations 

In this section we derive the \J~D equations by approximating the kinetic energy 
of vertical motion in Hamilton's principle for a multilayer ideal fluid. Such a system 
consists of N homogeneous fluid layers of densities pi,i e [1, . . . N] , where i = 1 is the 
top layer and i = N is the bottom layer. Thus, for stable stratification, p i+ \ > p t . The 
i-th layer has a horizontal velocity, Ui and thickness, A; the interface between the i- 
th and i— 1-th layer is at depth hi = —b+J2f=i A, for a prescribed bathymetry, b(x, y). 
We assume columnar motion within each layer (horizontal velocity independent of 
vertical coordinate, z). Incompressibility then implies that vertical velocity is linear 
in z. Under this ansatz and after a vertical integration, the Lagrangian for Euler's 
fluid equations with a free surface becomes, 
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In this Lagrangian, the w 4 should be viewed as functions of layer velocities, layer 
thicknesses and their partial deri vatives. The fluid momentum equations are ob- 
tained from Euler-Poincare theory dHolm et al.W 19981) as 
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The system is closed by the layer continuity equations, 

^ + V ■ Dim = 0. 
at 

To obtain the MGN equations one sets Wi in 02.1ft equal to the vertical component of 
the fluid velocity, represented in terms of the material time derivative as 



(2.2) 



N 



Wi = (hi — z)V • Ui + Ui ■ Vhi — V ■ Dj 



(2.3) 
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In each layer, we define a vertical material coordinate, < s t < l,by,Sj = (h l — z)/D l . 
Then, by using ds i /dt l + wds l /dz = with d/dt l = (d/dt + u t ■ V), we note 



Wi = D, 



ds t 
dU 



(2.4) 



To obtain the new \f~D equations, we replace Wi in ( 12. ID with a different linear 
approximation of the vertical motions within each layer, as follows 
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This approximation introduces a set of N new length scales, di, the far field fluid 
thicknesses. The final quantity in ( 12.51 ) is a vertical fluid velocity, expressed in the 
convective representation (Hol m et a/J|l986l) . The spatial and convective represen- 
tations of fluid dynamics are the analogues, respectively, of the spatial and body 
representations of rigid body dynamics on 50(3). This analogy arises because the 
configuration spaces for fluid dynamics and for rigid bodies are both Lie groups. In 
both cases, spatial velocities are right-invariant vector fields, while the convective, 
or body, velocities are the corresponding left-invariant vector fields. 

On taking variations, the final equations arising from Hamilton's principle for 
the new Lagrangian may be written compactly in terms of a shallow water equa- 
tion, plus additional nonlinear dispersive terms that represent a non-hydrostatic 
pressure gradient, 



dui _ _ 

— + Ui ■ V« ; = -g\7 
at 



3=1 Pl 



V 



dj Fj v — -s 



(2.6) 



hydrostatic pressure non-hydrostatic pressure 



F, 



d di d 
dtD~dt 



[2hi + lu+i] , G, 



d di d 
dtDldt 



[hi + h l+ i] 



E, 



dhi 



dhi dh 



dhu 



dt J dt dt V dt 
In comparison, the MGN equations have a similar form, but the non-hydrostatic 
forces no longer correspond to the gradient of a pressure. Instead, the MGN equa- 
tions take the form 
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Having arisen from Hamilton's principle, the \JD and MGN equations may both 
be written in Hamiltonian form. In addition, since their Lagrangians are both in- 
variant under particle relabelling, the two sets of equations each possesses the cor- 
responding materially-conserved quantities, 



dq 

~dt 
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The quantity is the potential vorticity in the z-th layer. For the Vd equations, 



q, t = z ■ (V x Ui) /Di, 

just as for the unmodified shallow water equations. Hence, equation ( 12.61 ) admits po- 
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tential flow solutions for which V x Ui = 0. In contrast, the MGN potential vorticity 
contains higher derivatives of ui. 
When restricted to a single layer, the Lagrangian for the \f~D equations ( 12. 6D is 



D\u\ 



(D-bf 



dx Ay, 



and the single-layer VD motion equation becomes 

du _ _ , „ . , d 2 



dt 



u-Vu= -gV (D - b) 
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The single-layer Lagrangian con tains a term that coincid es with the Fisher-Rao 
metric in probability theory (See IBrodv & Hughstonl dl998l) for an excellent review 
of the subject.) This feature suggested the name a -/D equations" for the new system. 



3. Linear Dispersion Analysis 

This section shows that the \f~D equations in H2.6D possess the same linear disper- 
sion relation at rest as MGN, but that unlike MGN th ey remain linearly well-p osed 



sion relation at rest as IvHjtJN, but mat unlike ivilxlN tney remain linearly well-posec 
when background shear is present. As discussed in iLiska & Wendroff ( 1997b . lin 



ear well-posedness requires that the phase speed remain bounded as k — > oo for all 
background shear profiles. In what follows we specialise to the case N = 2, and im- 
pose a rigid lid constraint, hi = 0, through an additional barotropic pressure term, 
although the analysis generalizes to an arbitrary number of layers and to the free 
surface case. Linearizing equations J2.6D and ( ]2.2t > for far field depths, d\, d 2 and 
background velocities, Ui, U 2 produces the following dispersion relation 

Pid 2 {Ui - \f + p 2 di (U 2 - A) 2 + (A 2 /c 2 d 1( i2/3) (pidi + p 2 d 2 ) - (pa - pi)gdid 2 = 0. 

Here the phase speed is A = cj/k for frequency w and wave number k. Meanwhile 
under the same assumptions the equivalent CC equations linearize to 

Pid 2 (t/i - A) 2 (1 + d!fc 2 /3) + padi (U 2 - A) 2 (1 + d 2 k 2 /3) - ( P2 - Pl )gdxd 2 = 0. 

When the background state is at rest, U\ = U 2 = 0, both sets of equations exhibit 
the same dispersion relation, 

A _ ±\/{P2 - Pi)gdid 2 



\J P\d 2 + p 2 di + (fhd\ + p 2 d 2 ) did 2 k 2 

in which the phase speed of the high wave number modes converges towards zero. 
This is not surprising, since when linearized at rest the two measures of vertical 
motion coincide in ( 12. 3D and ( 12.51 *. with Wi = Wi. 

In the presence of background shear, U± ^ U 2 , the solution for the phase speed A 
of the rigid lid \f~D equations is 

-B ± y/5 2 - 4 (A + A 2 k 2 ) (Cp - (p 2 - p 1 )g~d^d 2 ) 
2(A + A 2 k 2 ) 
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A 



- (B + B 2 k 2 ) ± y/(B + B 2 k 2 ) - 4 (A + A 2 k 2 ) (C + C 2 k 2 - (p 2 - p 1 )gd 1 d 2 ) 

2(A + A 2 k 2 ) ' 

for k > 1. 



-B 2 ± y/B'j - 4A 2 C 2 
2A 2 



Here the coefficients A , B Q , Cp, etc. are functions of the mean depth ratio, di/d 2 



and of the shear (U 2 — U\). Liska & WendroflN 19971) showed that the discriminant 
B\ - A 2 C 2 is strictly negative in the CC case, thereby ensuring the existence of a root 
with growth rate, 3(Afc), which is positive and scales linearly with k. Thus, the CC 
equations are linearly ill-posed. In contrast, the magnitude of growth rates in the 
\f~D case are bounded from above, independently of wave number k. Thus, the \[D 
equations are linearly well-posed. Indeed, provided a Richardson-number condition 
is satisfied, that 



Ri := (p 2 - Pi) gd,d 2 {di + d 2 ) 2 I ( Pl dl + p 2 d\) \U 2 - U t \ 



> 1, 



(3.1) 



then the linearized system remains stable to perturbations at any wave number. 



4. Travelling Wave Solutions 

At this point, we have shown that the VD equations satisfy requirements (a) 
and (b) in the Introduction. We will now show that the Vd and MGN equations 
also satisfy requirement (c); that is, they share the same travelling wave solutions 
for any number of layers, in the absence of imposed background shear. We take a 
travelling wave ansatz for a choice of wave speed, c, 

Di(x,t) =: Di(X), Ui(x,t) =: Ui(X), with X := x — ct . 

Integrating the continuity equation ( 12.21 ) in a frame moving at the mean barotropic 
velocity implies for both the MGN and \[T) systems that 

u i = c(l-d i /D i ). (4.1) 

For such solutions we see from equations 02.3H and 02.5D that Wi and Wi coincide, 

~(z - hi) dDi dhi 



( h . - Z )^i + ( u . - c )^!i = -^1 



Di dX dX 



= Wi. (4.2) 



To obtain travelling wave solutions, one varies the Lagrangian subject to the rela- 
tion d4.ll >. Using relation ( 14. 2D between the vertical velocities implies that the travel- 
ling wave Lagrangians for the \AD and MGN equations coincide. Consequently, the 
travelling wave solutions of the two equation sets must also coincide. This finishes 
the derivation of the new equation set and the verification of the three require- 
ments set out in the Introduction. In the presence of background shear, however, 
(w t — Wi) ^ 0, and the coincidence of the travelling waves of the two models no 
longer holds. 



5. Numerical experiments (Single layer) 

We now present some results of numerical computation with the \[D equations 
for a single layer. The equations are discretized using standard finite volume tech- 
niques, without any filtering. In this section we compare results for the one-layer 
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FIGURE 1. Single-layer \/T5 and GN thickness profiles arising from lock-release runs. The ar- 
rows indicate directions of crest propagation. The leading wave profiles are essentially iden- 
tical for both the \[~D and GN equations. The smaller crests are also very similar. (Differences 
are emphasised by using a logarithmic scale.) 

\f~D and GN equations. These results show that the two systems have the same 
qualitative behaviour. 



We first consider a lock release experiment with periodic boundary conditions and 
physical parameters and initial conditions given by 



D(x, 0) = 1 + 0.5 [tanh (x - 4) - tanh (x + 4)] , u(x,0) = Q, x £ [—100, 100]. 

We integrated this state forwards in time for both the GN and \fl) systems. Figure 
Q] shows a snapshot of a train of travelling waves emerging from the lock release in 
each equation set. In both cases, the first wave is the largest and the wave profiles of 
the two systems track each other closely. This may have been expected, because the 
two systems share the same travelling wave solutions and possess identical linear 
dispersion relations for a quiescent background. 

5.2. Single layer interaction of solitary waves 

We now consider the interaction between pairs of solitary travelling wave solutions 
in the two systems. The V~D and GN systems have the same travelling wave so- 
lutions, but their PDEs differ; so one may expect to see differences in their wave- 
interaction properties. For both systems the travelling wave takes a sech 2 form, 



and the velocity profile is given by equation 04.1D . Figure [2] shows a snapshot of 
the results for the case of two equal amplitude, oppositely directed travelling waves 
after suffering a symmetric head-on collision. Both the \[D and GN descriptions 
show a near-elastic collision followed by small-amplitude radiation. The travelling 
wave behaviour is essentially identical. The main difference lies in the details of the 
slower linear waves that are radiated and left behind during the nonlinear wave 
interactions. Similar results were obtained in the case of asymmetric collisions. 



5.1. Single layer lock release 



d=l, .9 = 1, 
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FIGURE 2. Single-layer \/T5 and GN results are shown for symmetric head-on collision exper- 
iments. Layer thickness profiles are shown after the collision. The fully nonlinear collision is 
near elastic for both \fT) and GN, with near conservation of the initial amplitude, shown by 
the dotted line. The travelling wave behaviour is essentially the same, except for a very small 
phase shift, while the slow, small-amplitude radiation waves show a few minor differences. 

6. Numerical Experiments (Two Layers, Rigid Lid) 

We now consider numerical solutions of the t wo-layer rig id lid Vd equations in 
one spatial dimension. Unlike the result of Jo & Ch oi (2002) for the CC equations, 
there is no value of resolution at which the code for numerically integrating the \J~D 
equations becomes unstable. This might have been hoped, because the \[T) equa- 
tions are linearly well-posed. 



We next present the result of the equivalent lock-release experiment to that of Sec- 
tion 15.11 for the two-layer rigid lid equations in lock release configuration. For this, 
we choose d± = 0.1, d 2 = 0.9, g = 1, and 

l>i(x,0) = di +0.3 [tanh(x- 4) - tanh (x + 4)] , D 2 = d x + d 2 - D lt x G [-100,100]. 

The initial configuration is taken as a depression of the fluid interface, because CC 
or \f~D travelling waves extend from the thinner layer into the thicker one. Since 
our motivating interest was in modelling waves on the ocean thermocline, we will 
concentrate here on the situation in which the top layer is thinner than the bottom 
layer. We show results for the following velocity profiles: 

(a) No shear: ui(x,0) =0,n 2 (x,0) = 0, and x G [-100,100]. 

(b) Background shear flow: u\(x, 0) = 0.01, u% = -DiMi/_D 2 ,amd x G [-100,100]. 
From A3.1H the Richardson number for the case with shear is Ri = 5.02, and the flow 
is still linearly stable for all wave numbers. 

The result of the integration is shown in Figure[3] Unlike the single-layer case, the 
linear long wave speed in the two-layer case is close to the wave speed of the result- 
ing travelling waves. In Figure |4] we compare the leading short wave disturbances 
of one lobe in the shear- free lock-release experiment to the calculated solitary trav- 
elling wave solutions in the same parameter regime, as denned by the maximum 
amplitude of the disturbance. For the three largest waves the differences are indis- 
tinguishable, while the smaller waves, which are still interacting with the long wave 
signal are wider than the equivalent true travelling wave. A similar result holds for 
the asymmetric crests generated in the presence of shear and for the respective ve- 



6.1. Two layer rigid lid lock release experiments 




FIGURE 3 . Two-layer rigid lid lock release experiments are compared, with and without back- 
ground shear of Richardson number Ri = 5.02. The fluid interface is plotted at times t = 0, 
600, 1200. Here 3 = 1, pi/p 2 = 0.995, di/hi = 0.1. The fluid was initialised with a sym- 
metric tanh profile in thickness. Filled arrows denote the directions of travel of the waves, 
while unfilled arrows indicate the background shear. The run without shear (upper figure) 
generated a left-right symmetric distribution of wave crests. In contrast, the run with back- 
ground shear (lower figure) generated an asymmetric distribution of wave peaks. The lower 
figure shows lock release for an upper layer moving with initially constant rightward velocity. 
Consequently, there was a net transfer of momentum toward the right, which significantly 
affected the distribution of wave crests. 



locity profiles. The leading order solution is thus well approximated as a train of 
independent solitary travelling waves propagating in order of amplitude. 

6.2. Interaction of two layer solitary waves with rigid lid 

Finally we consider wave interaction experiments in the two-layer case. Here we 
initialize with two numerical solutions of the travelling wave problem with oppo- 
sitely directed velocities. Results for these two e xperiment s are similar to those in 
the single layer case and also similar to those in I Jo & Ch oi (2002) in the case with- 
out shear. Figure [5] shows our numerical solution following the wave collision in the 
case of equal magnitude wave speed with and without a background shear. We again 
see nearly elastic collisions, after which the waves approximately regain both their 
shape and amplitude. In the case with shear, the wave profiles are no longer sym- 
metric. Defining upstream and downstream direction in terms of the background 
flow in the thicker layer, we see that the upstream gains amplitude and narrows 
compared with the shear-free wave profile at the same wave speed, while the wave 
travelling upstream loses amplitude and widens. Nevertheless, the waves are seen 
to nearly conserve amplitude, and hence momentum through the collision. 
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FIGURE 4. Close up of two-layer wave lock-release experiment in Figured] The corresponding 
travelling wave profile (solid line, as defined by the maximum wave amplitude) is plotted over 
the numerical lock-release solution at t=1300 (dashed line). In the three leading waves, the 
differences from the corresponding solitary wave profiles are indiscernible. 
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FIGURE 5. The figure shows two-layer head-on collision experiments with a rigid lid. In the 
top figure there is no background shear, while in the bottom figure the Richardson number is 
Ri — 502. The interface depth is shown after the symmetric head-on collision of two travelling 
waves of the same phase speed. Physi cal p arameters are, g = 1, pi/p-i = 0.995, di/fri = 0.1. In 
both figures the phase velocity is c/^/ghi = ±0.027 and the initial amplitudes are indicated 
by dotted lines. Filled arrows denote the directions in which the solitary waves are travelling, 
while unfilled arrows indicate the background shear. In both cases, the collision was found 
to be nearly elastic, and the final conditions were essentially mirror reflections of the initial 
conditions, with only slight changes in the amplitudes. 



7. Summary 

We have introduced a system of Hamiltonian multilayer water wave equations 
that is linearly well-posed in the presence of background shear but possesses the 
same travelling wave solutions as the MGN equations, whose ill-posedness had pre- 
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viously caused difficulties in their numerical integration. The new system also has 
been found numerically to generate fast-moving trains of large-amplitude coherent 
waves that exhibit ballistic, nearly-elastic, nonlinear scattering behaviour amidst a 
background of slow, small-amplitude, weakly interacting or linear wave radiation. 
Future work will explore wave-wave and wave-topography interactions in the pres- 
ence of vertical shear between the layers, as well as two-dimensional multilayer 
wave interactions and wave generation by flow over topography by using the new 
well-posed V~D equation set. 

We are grateful to W. Choi, F. Dias and R. Grimshaw for fruitful discussions on 
this topic. J. R. Percival was supported by the ONR grant N000 14-05- 1-0703. 
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